function  [Error_i, CS_i, W_i, W_experience_i,  Error_avg, CS_avg, W_avg, W_experience_avg, ext_i, ext_avg, ext_experience_i, ext_experience_avg, P_decision_j, P_experience_j,gvt_i,gvt_avg,SW_i,SW_avg,gvt_experience_i,gvt_experience_avg,SW_experience_i,SW_experience_avg]=Welfare_k_clogit_EEgap(param,param_FE,param_FEdemo,id_ij,J,p_elec,t_pigou,d_ext,elec,lifetime,r,MCPF)

% Welfare formula:
% Delta CS = 1/eta *([log (sum_j U_j) +sum_j P_j(U^e-U)]_post_policy - [log (sum_j U_j) +sum_j P_j(U^e-U)]_pre_policy)

global state_i week_i T S MT MTvector K WB...
price price_denom rebate_estar_denom estar_denom kwh prod_denom...
id_i id_j id_m id_n  demo pidnest choicenest_ij brandweek_denom ...
type_id2_3_demo1_num type_id2_3_demo1_denom type_id3_demo1_num type_id3_demo1_denom AV_demo1_num AV_demo1_denom  ice_sc_demo1_num ice_sc_demo1_denom...
type_id2_3_demo2_num type_id2_3_demo2_denom type_id3_demo2_num type_id3_demo2_denom AV_demo2_num AV_demo2_denom  ice_sc_demo2_num ice_sc_demo2_denom...
type_id2_3_demo3_num type_id2_3_demo3_denom type_id3_demo3_num type_id3_demo3_denom AV_demo3_num AV_demo3_denom  ice_sc_demo3_num ice_sc_demo3_denom...
type_id2_3_demo4_num type_id2_3_demo4_denom type_id3_demo4_num type_id3_demo4_denom AV_demo4_num AV_demo4_denom  ice_sc_demo4_num ice_sc_demo4_denom...
type_id2_3_demo5_num type_id2_3_demo5_denom type_id3_demo5_num type_id3_demo5_denom AV_demo5_num AV_demo5_denom  ice_sc_demo5_num ice_sc_demo5_denom...
type_id2_3_demo6_num type_id2_3_demo6_denom type_id3_demo6_num type_id3_demo6_denom AV_demo6_num AV_demo6_denom  ice_sc_demo6_num ice_sc_demo6_denom;


 %Parameters
  alpha=[0;param_FE(1:J-1,1)];
  eta=param(1);
  eta_experience=param(2);
  theta=param(3);
  tau=param(4);
  pi=param(5);
  beta=param_FEdemo;
  
  FEprod_denom=alpha(prod_denom); 
 
 % kwh  
 kwh_denom=kwh(id_ij);
  
 % elec
    p_elec_ext=elec+p_elec+t_pigou;
  	elec_J=repmat(p_elec_ext,1,J);
    elec_J=elec_J';
    elec_denom=elec_J(id_ij).*kwh_denom;
    
       
 %Computations  
 
 U_decision= FEprod_denom+eta*price_denom +theta*elec_denom...
 				+tau*estar_denom + pi*rebate_estar_denom... 				
 				+beta(1)*type_id2_3_demo1_denom +beta(2)*type_id2_3_demo2_denom + beta(3)*type_id2_3_demo3_denom...
  				+beta(4)*type_id2_3_demo4_denom +beta(5)*type_id2_3_demo5_denom + beta(6)*type_id2_3_demo6_denom...
				+beta(7)*AV_demo1_denom +beta(8)*AV_demo2_denom + beta(9)*AV_demo3_denom...
  				+beta(10)*AV_demo4_denom +beta(11)*AV_demo5_denom + beta(12)*AV_demo6_denom;
 
 expU_decision=exp(U_decision);
  				
 discount_fac=1/(1+r); 
 LLC=discount_fac*(1-discount_fac^lifetime)/(1-discount_fac); 
 eta_LLC=eta_experience*LLC; 				
 
 U_experience= FEprod_denom+eta_experience*price_denom +eta_LLC*elec_denom...
 				+tau*estar_denom + pi*rebate_estar_denom... 				
 				+beta(1)*type_id2_3_demo1_denom +beta(2)*type_id2_3_demo2_denom + beta(3)*type_id2_3_demo3_denom...
  				+beta(4)*type_id2_3_demo4_denom +beta(5)*type_id2_3_demo5_denom + beta(6)*type_id2_3_demo6_denom...
				+beta(7)*AV_demo1_denom +beta(8)*AV_demo2_denom + beta(9)*AV_demo3_denom...
  				+beta(10)*AV_demo4_denom +beta(11)*AV_demo5_denom + beta(12)*AV_demo6_denom;				
  					
  expU_experience=exp(U_experience);
    
  U_decision_ij = sparse(id_i,id_j,U_decision,id_m,id_n)';
  num_decision = sparse(id_i,id_j,expU_decision,id_m,id_n)';
  denom_decision=sum(num_decision,2);

  U_experience_ij = sparse(id_i,id_j,U_experience,id_m,id_n)';
  num_experience = sparse(id_i,id_j,expU_experience,id_m,id_n)';
  denom_experience=sum(num_experience,2);
  
  P_decision=num_decision./repmat(denom_decision,1,J);
  P_experience=num_experience./repmat(denom_experience,1,J);
  
  P_decision_j=sum(P_decision,1)/MT;
  P_experience_j=sum(P_experience,1)/MT;
  
  Error_i=1000*abs(1/eta_experience)*sum(P_decision.*(U_experience_ij - U_decision_ij),2);
  CS_i=1000*abs(1/eta_experience)*(log(denom_decision));
  W_i=CS_i+Error_i;
  W_experience_i=1000*abs(1/eta_experience)*(log(denom_experience));
  
  Error_avg=sum(Error_i)/MT;
  CS_avg=sum(CS_i)/MT;
  W_avg=sum(W_i)/MT;
  W_experience_avg=sum(W_experience_i)/MT;
   
  ext_ij=LLC*kwh*d_ext*1000;
  ext_i=sum(P_decision.*ext_ij',2);
  ext_avg=sum(ext_i)/MT;
  ext_experience_i=sum(P_experience.*ext_ij',2);
  ext_experience_avg=sum(ext_experience_i)/MT;
  
  gvt_ij=LLC*kwh*t_pigou*1000;
  gvt_i=sum(P_decision.*gvt_ij',2);
  gvt_avg=sum(gvt_i)/MT;
  gvt_experience_i=sum(P_experience.*gvt_ij',2);
  gvt_experience_avg=sum(gvt_experience_i)/MT;
  
  SW_i=W_i-ext_i+MCPF*gvt_i;
  SW_avg=sum(SW_i)/MT;
  
  SW_experience_i=W_experience_i-ext_experience_i+MCPF*gvt_experience_i;
  SW_experience_avg=sum(SW_experience_i)/MT;
end 











